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Abstract. We implement and further refine the recently proposed method (Kashlinsky, Hernandez-Monteagudo & Atrio- 
Barandela, 2001 - KHA) for a time efficient extraction of the power spectrum from future cosmic microwave background 
(CMB) maps. The method is based on the clustering properties of peaks and troughs of the Gaussian CMB sky. The procedure 
takes only ^[f(v)] 2 N 2 steps where f(v) is the fraction of pixels with |<5r| > v standard deviations in the map of N pixels. 
We use the new statistic introduced in KHA, f v , which characterizes spatial clustering of the CMB sky peaks of progressively 
increasing thresholds. The tiny fraction of the remaining pixels (peaks and troughs) contains the required information on the 
CMB power spectrum of the entire map. The threshold v is the only parameter that determines the accuracy of the final spec- 
trum. We performed detailed numerical simulations for parameters of the two-year WMAP and Planck CMB sky data including 
cosmological signal, inhomogeneous noise and foreground residuals. In all cases we find that the method can recover the power 
spectrum out to the Nyquist scale of the experiment channel. We discuss how the error bars scale with v allowing to decide 
between accuracy and speed. The method can determine with significant accuracy the CMB power spectrum from the upcoming 
CMB maps in only ~ (10~ 5 - 10~ 3 ) x N 2 operations. 
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k*" ■ 1 . Introduction. power spectrum, on sub-degree scales the interaction between 
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Oj ■ were inside the horizon at the last scattering epoch. The CMB on a final set of the cosmo i ogica i parameters (Q total , Q A , Q b , h) 

fluctuations on these scales carry a signature of causal pro- an d also serve to validate the cosmological CDM paradigm (see 

cesses during the last scattering and thereby provide a very Hu & Dode i son 2002 for a recent review). It is then important 

important constraint on the physics of the early Universe and tQ meafjure the sub . degree structure G f the CMB with high ac- 

the models for structure formation. The most popular of these _,„„_,. 

A A curacy. 

models is the cold dark matter (CDM) set of models based 
on the inflationary paradigm for the evolution of the early 

Universe. The models are very appealing, not only because of After its first year of operation, the WMAP experiment 

their relative simplicity, but also because they provide a clear- has measured the Cosmic Microwave Background (hereafter 

cut set of predictions that can be verified by observations. One CMB) at five different frequencies with the maximal angu- 

(and perhaps the most critical) of these predictions is the sub- lar resolution of ~ 0.21° and sensitivity close to 175 yuK per 

degree structure of the CMB anisotropics. In the framework of 7' pixel, (Bennet et al. 2003a). With these sensitivity and an- 

CDM models with adiabatic density fluctuations the structure gular resolution levels, the cosmological parameters as Q , 

of the CMB power spectrum reflects the linear physics of sound &a, ^baryon, H , n, T reio were determined with high accuracy 

waves and initial density perturbations. While on scales outside (Spergel et al 2003). These measurements should be further 

the horizon at de-coupling the CMB field preserves the initial improved with the planned ESA's mission Planck, where maps 

will contain up to 10 7 pixels and extend to higher angular res- 

* Formerly at: Facultad de Ciencias, Universidad de Salamanca, Pza olution. Angular resolution and noise are both important for 

de la Merced, s/n, 37008, Spain determining the cosmological parameters from the CMB. 



the photon fluid and matter leads to a series of acoustic peaks. 
The sub-degree structure of the CMB probes linear scales that The rdative heigm% width and spacing of these peaks depend 
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A major challenge to understanding current and future 
CMB measurements is to find an efficient algorithm that can 
reduce these enormous datasets: N ^ 10 5 pixels in balloon ex- 
periments, =a 3 x 10 6 for WMAP, to 5 x 10 7 for the Planck 
HFI data. (For comparison the COBE DMR data analysis was 
based on only 4,144 pixels). Traditional methods require in- 
verting the covariance matrix and need ~ A^ 3 operations making 
them unfeasible for the current generation of computers. Thus 
alternatives have been developed for estimating the CMB mul- 
tipoles from Gaussian sky maps: Tegmark (1997) proposed a 
least variance method that yields C/'s directly from the temper- 
ature map in 0(N 2 ) operations, while Oh, Spergel & Hinshaw 
(1999, hereafter OSH) perform an iterative maximization of 
the likelihood of the temperature map, also in 0(N 2 ) opera- 
tions. The first year WMAP results were analyzed with the 
OSH method (Hinshaw et al 2003). Bond, Jaffe & Knox (2000) 
and Wandelt, Hivon & Gorski (2001) concentrate on the statis- 
tics of C/'s, once the temperature map has been Fourier trans- 
formed. Hivon, et al. (2002) studied the likelihood of the power 
spectrum as obtained by direct Fast Fourier Transform of the 
available portion of the sky. Though it requires 0(N 3 ^ 2 In N) 
operations, the accuracy of their results depends on the fraction 
and geometry of the sky covered. A different approach consists 
in computing the correlation function directly from the data in 
0(N 2 ) operations. Smoot et al. (1992) used this type of analy- 
sis for the first year release of the COBE data. More recently, 
Szapudi et al. (2001), Szapudi, Prunet & Colombi (2001) have 
developed it for mega-pixel CMB data sets. The last contribu- 
tion in this field comes from Penn (2003), who applies existing 
(Padmanabhan, Seljak & Penn 2002) iterative algorithms on 
CMB data analysis, and manages to estimate to power spectra 
in 0(N log N) operations. 

A different method to compute the CMB power spectrum 
in a fast and accurate manner was proposed by us (Kashlinsky, 
Hernandez-Monteagudo & Atrio-Barandela, 2001, hereafter 
KHA). The method exploits Gaussian properties of the CMB 
and noise fields and uses high peaks (and troughs) of the CMB 
field whose abundance is much smaller than the total number of 
pixels, N, and whose correlation properties are strongly ampli- 
fied in a way that depends on the underlying power spectrum. 
This simultaneously achieves two important goals: reducing 
the number of computational steps for analysis and the good 
accuracy of the measured parameters. KHA have shown that 
the tiny fraction of the remaining pixels (peaks and troughs) 
contains the required information on the CMB power spectrum 
in the small scales. The peaks also trace, by default, the pix- 
els with high signal-to-noise ratio and keep most of the infor- 
mation about the power spectrum of the signal. Although this 
method may not provide as accurate results as other more time 
consumings methods, it provides a new, non-standard and in- 
dependent tool to unveil the sub-degree structure of the CMB. 

In this paper we present detailed numerical simulations 
with the application of the KHA method to WMAP and Planck 
datasets in the presence of realistic components, such as inho- 
mogeneous noise and non-Gaussian features expected from the 
Galactic foreground residuals. We show that with this method 
we can determine the CMB power spectrum outside the beam 
(Z * 640 for WMAP and higher for Planck) in only < 10~ 3 N 2 



operations. We show that for the projected two-year WMAP 
noise levels our error bars at each I are, on average, compara- 
ble to OSH, but have larger correlations which may be a small 
price to pay for the significantly shorter computational time. 
For both WMAP and Planck parameters the KHA method is 
also faster than the direct computation of the CMB power spec- 
trum in 0(N 3/2 log N) ~ N L61 steps. 

The structure of this paper is as follows: In Sect. 2 we 
briefly review the KHA formalism and in Sect. 3 we discuss its 
numerical implementation. Sect. 4 deals with application of the 
method to both idealized and realistic WMAP data. We show 
there that the method is immune to noise inhomogeneities and 
non-Gaussian features of Galactic foregrounds. Sect. 5 follows 
with results for application of the KHA method to Planck and 
in Sect. 6 we present our main conclusions. 

2. Mathematical formalism: an overview 

For completeness we give a brief overview of the KHA method; 
more details are in KHA. 

The CMB sky is expected to be highly Gaussian (Komatsu 
et al 2003) and this property, Eq. ([0 below, is also widely 
used in standard maximum likelihood methods. For a Gaussian 
ensemble of N data points (e.g. pixels) describing the CMB 
data 6 = T - (T) one will find a fraction f(v) =erfc(v/ y2) 
with \6\ > vcr, where cr 2 - (6 2 ) is the variance of the field 
and erfc is the complementary error function. The fraction of 
peaks, /(v), is a rapidly decreasing function for vl 1 and is e.g. 
/(v) = (4.5, 1, 0.1) x 10~ 2 for v = (2, 2.5, 3) respectively. The 
joint probability density of finding two pixels within d6\^ 2 of 
(5i,2 and separated by the angular distance 8 is given by the bi- 
variate Gaussian: 



1 



exp(-itf-C _1 -6) 



(1) 



where C is the covariance matrix of the temperature field. We 
model the temperature correlation function (or matrix, if de- 
fined on a set of pixels) as 



C(0 y ) = C CMB (%) + (NiNj) 



(2) 



where C CMB (0 l7 ) and (NiNj) are the CMB and the noise cor- 
relation matrices, respectively. We further define the total vari- 
ance of the map as Co = C CMB (0) + (A/ 2 ). The case of inhomo- 
geneous noise, i.e., noise whose variance varies across the sky, 
will be discussed in Sect. 4.2. Note that the power spectrum 
of the map is nothing but the Legendre transform of C{6), (see 
Eq. (0 below). 

The distribution of peaks of a Gaussian field is strongly 
clustered (Rice 1954, Kaiser 1984, Jensen & Szalay 1986, 
Bardeen et al. 1986, Kashlinsky 1987). Their angular cluster- 
ing can be characterized by the 2-point correlation function, £ 
describing the excess probability of finding two events at the 
given separation. The correlation function of such regions is: 

{ v (8) = ^ Lj ^ ts 1=A V (— )(3) 

[2 m p{6)d5f Q 
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where A v is evaluated in detail in KHA to be: 



A V (X) 



V2 y k=l 



T^-t(^) 



(4) 



Here H n {x) = (-)" exp(x 2 )d n / dx" exp(-jc 2 ) is the Hermite 
polynomial; H-i(x) = exp(x 2 )erfc(x). 

At each angular scale the value of £ v for every v is deter- 
mined uniquely by C at the same angular scale. In the limit 
of the entire map (v=0) our statistic is £ v =0 and our method 
becomes meaningless; the new statistic has meaning only for 
sufficiently high v. One should distinguish between the 2-point 
correlation function, we directly determine from the maps, 
and the commonly used statistics in CMB studies, the temper- 
ature correlation function, C. 

KHA thus suggested the following procedure to determine 
the power spectrum of CMB in only ^ f 2 (v)N 2 <K A^ 2 opera- 
tions: 

1. Determine the variance of the CMB temperature, Co, 
from the data in N operations; 

2. Choose sufficiently high value of v when /(v) is small 
but at the same time enough pixels are left in the map for robust 
measurement of £ v (#); 

3. Determine % v (9) in [f(v)] 2 N 2 operations. 

4. Finally, solve the equation A v (C/Co) = £ y (0) to obtain 
C CMB {9) and from it C/. This last step assumes that the exper- 
imental noise correlation function (NiNj) has been properly 
determined, at least in the set of pixels selected in the analysis, 
so that it can be subtracted from the total correlation function 
C(0y). 

Formally speaking this procedure would require 
j[f(v)] 2 N 2 operations because of the symmetry in count- 
ing pairs. For simplicity we will be referring to this as 0(f 2 N 2 ) 
method implying a gain factor of [/(v)] 2 compared to other 
0(N 2 ) methods. 

3. Numerical implementation 

3.1. Modeling the CMB sky 

Maps of the CMB sky were generated using the hierarchical 
equal area isolatitude pixelization of the celestial sphere imple- 
mented in HEALPix 1 . Pixels have equal area and are arranged 
in "constant latitude rings". Maps can be constructed with vary- 
ing resolution, the number of pixels given by N = 12xAf 2 de , be- 
ing A^side the number of times in which each side of a pixel will 
be divided in two, starting from a given initial configuration. 
For a value of N s id e of 512, one generates a map of 3,145,728 
pixels of size of seven arcminutes. In the case of PLANCK, we 
used A^side = 1024 or 12,582,912 pixels of 3.5 arcminutes size. 
The CMB spectrum was obtained using the CMBFAST code 
(Seljak & Zaldarriaga, 1996). 

For WMAP and PLANCK, we performed two types of 
simulations including: (1) cosmological signal plus white 
noise, and (2) cosmological signal, foregrounds and inho- 
mogeneous white noise. We chose three different thresholds: 
v = 2.0,2.27,2.525 for WMAP and v = 2.55,2.78,3.0 for 



Planck. At each threshold, the pixels selected are (on average) 
f(v)N and their number doubles with decreasing threshold. The 
thresholds for WMAP and Planck were chosen such that they 
select the same number of pixels in both experiments. In this 
way, we can analyze the effect of pixel number on the accuracy 
of the method, i.e., we can estimate how the multipole error 
bars scale with threshold. 

From the peak spatial distribution, % v (9) was computed on 
a grid of 31,415 equally spaced bins. The power spectrum is 
computed using a Gauss-Legendre integration: 



Ci = 2n 



f 

Jo 



d9sm9C(9)Pi(cos9). 



(5) 



Accurate integration requires the evaluation of C( ff) at the roots 
of the Legendre polynomial out to a maximum order l mdx (Press 
et al 1992) and this property was used by Szapudi et al (2001), 
Szapudi, Prunet, & Colombi (2001) and KHA. The value of 
^max = 800 was chosen as a compromise between the beam 
scale and the pixel angular scale, in order to prevent the first 
root of the Legendre polynomial from being within the pixel 
size where there is no information on % v (9). Prior to inverting 
Eq. % v (ff) was smoothed using a Gaussian of width 4' cen- 
tered on the roots of the Legendre polynomials. We verified 
that other average schemes lead to essentially identical results. 

Fig. CJi) shows £ v=2 (9) for a SCDM model, with fi tota i = 1, 
Q A = 0, O b = 0.04, H = 55 km s -1 Mpc and scalar per- 
turbation spectral index n$ = 1 . The solid line corresponds to 
the theoretical prediction and the dashed line is a smooth av- 
erage obtained from the data. The agreement is very good out 
to 9 ~ 10°, allowing an accurate reconstruction of the power 
spectrum in the range of interest: 30. 

3.2. Accuracy of £, and correlations 

There are numerous ways to estimate % v {9) from the data, each 
way coming with its own uncertainty and systematics. Ideally, 
the uncertainty in % v {9) should reach 1 / -sjN pa i rs but in practice 
systematic and other effects become dominant (Peebles 1980). 
Because the KHA method involves a complex and non-linear 
algorithm for inverting C/'s from % Y {9), the accuracy of the re- 
sultant CMB power spectrum may depend sensitively on the 
details, uncertainties and correlations of a particular estimator 
aff v . 

The ultimate accuracy with which the CMB power spec- 
trum can be determined at each / is given by the so called op- 
timal variance containing both the cosmic variance and the in- 
strument noise (Knox 1995): 



cov(0 = 



2 An{N 2 )BT 2 

-i/ [Ci + — — — ]. 

^(2Z+l)/ sky L N 



(6) 



In this equation, N is the number of pixels, f^ y is the frac- 
tion of the sky covered by the experiment and Bi is the window 
function due to the finite beam resolution. This uncertainty also 
limits the accuracy with which one can determine the correla- 
tion functions C and The latter are given by: 



HEALPix's URL site: {http://www.eso.org/science/healpix/) 
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Fig. 1. (a) g v=2 for the SCDM model with a Harrison- 
Zel'dovich power spectrum, obtained from a single realiza- 
tion of WMAP 94GHz channel. Solid line: theoretical estimate 
(Eq. 0)- Dots: raw data of £ v =2 evaluated at 31,415 angular 
bins distributed uniformly from 9 — 0° to 180°. Dashed line: 
result of filtering the raw data with a Gaussian of 4' width at the 
roots of the Legendre polynomial of order 800. (b) The shaded 
area represents the l<x optimal variance error bar (eq (6)) for 
the 94 GHz WMAP channel. Dashed and solid lines have the 
same meaning as in (a). 



The sampling variance associated with the low number of 
pairs - compared with methods based on the correlation func- 
tion - is much smaller that the cosmic variance on £ v . Fig. 
shows -^V pa j rs O£ov which is the ratio of the optimal variance 
uncertainty on ^ to the statistical uncertainty, 1 / JNp^. The 
lines correspond to various y=1.5 (dashes), 2 (dots) and 3 (solid 
line). The ratio was computed for WMAP number of pixels. 
The plot shows that for any mega-pixel CMB map the method 
can determine the new statistic, £„, out to the angles of interest 
optimally even for threshold as low as v = 1.5. This is the rea- 
son why our results are not sensitive to the particular estimator 
used in calculating £ v . 



100 



The shaded region in Fig. Q shows an example of the optimal 
variance uncertainty in £ v= 2(#). (Note that it is a function of v). 

In KHA, we evaluated £ v as the ratio of number of peak 
pairs at a given angular scale with respect to a poissonian cata- 
log: 

$ Y {ff) = DD/RR-l. (8) 

At each angular separation, DD is the number of pairs of peaks 
and RR the number of pairs of f(v)N points randomly located 
on the sphere Landy & Szalay (1993) argue that this estimator 
is neither optimal nor unbiased. They proposed a more accurate 
estimator defined by: 

£,(0) = (DD + RR- 2DR)/RR, (9) 

with DR is the number of pairs given a cross-correlation be- 
tween the f(v)N peaks of the CMB map and the same number 
of points Poisson distributed on the sky. We performed sim- 
ulations for WMAP as described in Sect. 4 and computed £ v 
using both estimators. We found that, for the range of interest 
(0 < 10°), both estimators show scatter equally close to Eq. Q. 
Hence, in what follows we shall use Eq. to be consistent 
with KHA. 

In Fig. {n>) we plot the 1-cr errors in £ v= 2 for the 94 GHz 
WMAP channel given by the noise and cosmic variances. The 
smooth average obtained from the data (dashed line) is in good 
agreement with theoretical value. Therefore, we will expect 
that the radiation power spectrum obtained from ^ v= 2 as de- 
scribed below will be very close to that of the input model. 



0.1 1.0 10.0 

0, (degrees) 

Fig. 2. The ratio of the optimal variance uncertainty on £ v to the 
statistical uncertainty, 1 / JN^s plotted vs the separation angle 
for the resolution of the WMAP 94GHz channel. The solid line 
corresponds to v = 3, dotted to v = 2 and dashes to v = 1.5. 



At 01 10° - 20° the value of £,(0) is small and its value 
becomes dominated by shot noise. Hence, we restrict the anal- 
ysis to < 10° - 20° by introducing a taper function that cuts 
out the contribution of the correlation function for 01 10°. As 
a result (a) we will not recover very accurately multipoles be- 
low / a 30 and (b) tapering will introduce additional correla- 
tions among the different C/'s. The first limitation is irrelevant, 
since one can always degrade the map to smaller resolution 
and apply standard techniques to recover the power spectrum 
at I < 20 - 40. The second limitation is common to all methods 
that compute the power spectrum by means of the correlation 
function or in the presence of Galactic (and other) cut. The in- 
trinsic multipoles need to be de-convolved from the tapering 
function, i.e., if Cj ntnnslc are the multipoles of the sky radiation 
power spectrum, and C\ are obtained using Eq. (0, then 

C, = J] C~ c T,-,> (10) 
where Ti is the Legendre transform of the tapering function F. 
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Fig. Q shows the Legendre transform of the tapering func- 
tion as function of Al = I — I' for Gaussian and top-hat tapering. 
The advantage of Gaussian tapering is quite obvious as it leads 
to significantly less prominent side-lobes. Hence, we adopted 
it in our computations. As the figure shows, it would lead to 
FWHM correlation width of only A/ ~ (a few) for tapering 
angles ; 10°. Correlations on the C/'s would be dominated by 
tapering and correlated errors on £ v . 



A Gaussian tapering 


A 







20 40 60 80 100 20 40 60 80 100 

M hi 

Fig. 3. Legendre transform of the tapering function (Eq. HOi 
plotted vs Al = I — I'. Left panel corresponds to Gaussian ta- 
pering, F(9) cc exp(-# 2 /6?) an d the right panel to the top hat 
tapering, F(9) = 1 for 9 < 9, and zero otherwise. Dashed, dot- 
ted and solid lines correspond to 9, = 18°, 12°, 8° respectively. 



tion and the power spectrum. Fig. @ summarizes the results 
for the threshold v = 2. In panel (a) the dotted line shows the 
power spectrum recovered from a single simulation, whereas 
the solid line gives the input power spectrum. Tapering angle 
was 9, = 12°. In (b), we show the average results for 500 re- 
alizations: the solid and dashed lines represent the input model 
and the average power spectrum of the 500 simulations, respec- 
tively. The shaded area shows the 1-cr error bar region for each 
C/ and the dashed lines show the optimal variance. Our method 
traces the acoustic peaks up to the beam scale, l\, elLm ~ 640. 
Its accuracy is roughly similar to the optimal variance cr ov for 
I < 300. The spectrum of Fig. has been obtained in 10~ 3 N 2 
operations, using a taper window of 12° FWHM. The accuracy 
can be improved by lowering the threshold v and hence increas- 
ing the computational time. 

To see the effect of tapering, in Fig. (0];) we plot the power 
spectra obtained for two different tapering angles for the same 
simulation as in Fig. ©). The solid line corresponds to 9, = 8° 
while the dashed line to 9, = 18°. Increasing the tapering angle 
leads to larger oscillations, but also decreases the correlations 
between different I. In Fig. 01) we plot the band power aver- 
age of the previous spectra: diamonds correspond to the taper- 
ing angle of 8° and triangles to 18°. The solid line is the input 
model and the dashed lines correspond to the optimal variance 
error bars as in (b). As expected, both results are almost identi- 
cal since all the information at high Fs is contained at angular 
scales smaller than 1°. 



4. Application to WMAP 

WMAP is observing at 5 frequency bands: 23, 33, 41, 61, and 
94 GHz, (Jarosik et al. 2003, Bennett et al. 2003a). The beam 
response at each band are given by a gain pattern, G, which 
is asymmetric and non-gaussian. From the solid angle beam, 
one can always define a FWHM, which, for increasing fre- 
quency channels, are equal to 0.82, 0.62, 0.49, 0.33 & 0.21 
degrees, (Page et al. 2003). The sky maps based on the full 
two-years of data are expected to have an rms noise of ^35 fiK 
per 0.3° x 0.3° pixel. By design, the noise will be essentially 
uncorrelated from pixel to pixel and it is expected to be highly 
Gaussian (Hinshaw 2000, Hinshaw et al. 2003). Due to the sky 
scanning strategy, the noise is reasonably uniform across the 
sky, except at the small regions near the ecliptic poles and at 
ecliptic latitude ~ 45° where the sensitivity will be somewhat 
higher. The WMAP radiometers produce raw temperature mea- 
surements that are the differences between two points on the 
sky separated ~ 140° . Since a given pixel i is observed with 
up to 1000 different pixels j, the covariance between any given 
pair of pixels (i, j) is much less than 1 % of the variance of pixel 
i. The noise covariance in the final sky maps will, by design, be 
very nearly diagonal. 

4.1. Homogeneous noise results 

First, we performed a set of 500 simulations for v = 
(2,2.27,2.525) including only cosmological signal and white 
noise. For each simulation, we computed the correlation func- 
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Fig. 4. Radiation power spectrum for v = 2 and different taper 
angles. In all panels the solid line shows the input power spec- 
trum, (a) Dotted line: power spectrum obtained from a single 
simulation with 9, = 12°. (b) Dotted line: power spectrum ob- 
tained from averaging 500 simulations, shaded area: l<x error 
region obtained from the same simulations. The dashed lines 
limit the l<x optimal variance error bars (Eq. (|6j). (c) The same 
as in (a) but for 9, = 8° (solid line) and 9, = 18° dashed line, (d) 
solid and dashed lines as in (b); the symbols displayed the band 
averaged power spectra with bandwidth AZ = 40 of the spectra 
plotted in (c): diamonds correspond to 9 t = 8° and triangles to 
9, = 18°. 
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The correlation matrix, Cy, among different multipoles is 
given by: 

etf . <jww_ i 

^(SCjXSC;,) 

with 6C1 = C\ — (Ci). The correlation coefficient matrix has 
been computed after using a Gaussian taper function F(6) = 
exp-(6>/6>,) 2 . It is highly diagonal with FWHM width Al ~ 10 
for 8, = 8°, independently of the value of the multipole /. At 
Al ~ 20 the correlation drops to the level of ~ 5%, indepen- 
dently of v. Outside the central diagonal strip the residual cor- 
relations are due to shot noise from the finite number of simu- 
lations. For consistency, we repeated the analysis with only 125 
simulations. The width of the diagonal remained the same, but 
the off-diagonal terms grew by a factor of two, consistent with 
Poisson statistics. 

For larger/smaller tapering angles, the correlation scale Al 
will be smaller/larger but, as demonstrated in Fig. (HJl), it 
will not change our estimate of the power spectrum. We also 
checked that our results did not change if we use the correla- 
tion function from the coarser (and computationally tractable) 
map from 6 = 10° out to 180° instead of tapering. The correla- 
tion on the final C/'s are dominated by the correlated errors in 
£ v at small angles. 

In Fig. 10 we compare the accuracy of our method for 
WMAP and v = 2 with that of OSH. We compute the ratio 
of the uncertainties in the multipoles recovered by both meth- 
ods. Our method gives a comparable precision to that of OSH 
but requires much fewer operations. On the other hand, it gives 
fewer independent data points. In direct methods, such as OSH, 
the necessary Galactic cut leads to C/'s that are correlated on 
scales of Al 2-3, while our method gives a correlation length 
of about Al ~ 10. 

In Fig. (|6j we show the variation of the error bar with band- 
width for the multipole at / = 200. Solid, dashed and dotted 
lines correspond to v = 2.5, 2.27 and 2, respectively. Due to the 
presence of this correlation among the multipole estimates, this 
cr C( ~ (AfT 1 ^ 2 behavior is obtained for All 10, the correlation 
scale introduced by our method. 

In Fig. <|7} we show cr C( at / = 200,300,400,500 for 
the three thresholds. Smaller error bars correspond to smaller 
thresholds, i.e., to larger number of peaks. A power-law fit 
using all multipoles gives a power law behavior of the form 
0"c, x N p . , with B = -0.405 ± 0.016. These relations can be 

L ' peaks' r 

used to estimate the amplitude of the error bar attached to each 
multipole for a wide range of values of v and Al. In particular, 
it can be used to find what values of v and Al are necessary to 
achieve a given degree of accuracy in the power spectra and the 
time required in the computation. 

We have shown that for a given experiment the accuracy of 
the C/'s computed spectrum depends mainly on one parameter: 
the peak threshold allowing to choose between larger accuracy 
or speed. 



2.0 




100 200 300 400 500 600 700 
Mutipole I 

Fig. 5. Comparison of the error bars of OSH and KHA methods 
for a CMB map as seen by the 94 GHz channel of WMAP. We 
plot the ratio (crq) KH A/(0"c,)osH- The uncertainties cr C( 's were 
computed for raw multipole estimates C/, without averaging. 



4.2. Including foregrounds and realistic noise 

We now generalize the method to include realistic foreground 
emission and inhomogeneous noise and will demonstrate that 
also in this case the KHA method gives similar estimates of the 
CMB power spectrum. 

4.2.1 . Foregrounds 

After the first year data release, the WMAP team made pub- 
lic foreground maps built after combining the five different 
frequency maps provided by the instrument, (Bennett et al. 
2003b). These foreground maps were constructed after apply- 
ing a Maximum Entropy Method using existing templates as 
priors, and the resulting model for galactic emission matches 
the observed emission to < 1%. Using this model for Galactic 
emission, we get that with the fiducial value of 20% for the 
amplitude of the foreground residuals above b — 10°, they are 
much smaller than the intrinsic CMB temperature fluctuations 
associated with peaks (ST 1 200^/K for v = 2). 

To model the effect of foregrounds, we used simu- 
lated maps provided to us earlier by G. Hinshaw of the 
WMAP science team (2002, private communication). These 
were produced by combining the Haslam 408 MHz map 
and the Schlegel, Finkbeiner & Davis (1988, hereafter SFD) 
1RAS/D1RBE 100 yum dust map. The Haslam map was used as 
the template for synchrotron emission and the SFD map as the 
template for dust and free-free. These maps were scaled to mi- 
crowave frequencies using the COBE DMR-based fits of these 
templates (with 7 degree resolution). The frequency by fre- 
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0.' 



0.0' 



I = 200 



will contribute negligibly to the angular power spectrum at 94 
GHz for I < 800. 




(90 GHz) fiK 



Fig. 6. Scaling of band power error bar <x C( with bandwidth Al. 
Solid, dashed and dotted lines correspond to v = 2.5, 2.27 and 
2, respectively. 



Fig. 7. Band power error bars as a function of the number of 
pixels used to compute the radiation power spectrum are shown 
for four values of I. Diamonds, triangles and squares corre- 
spond to Al = 30, 50, 70, respectively. 



quency fit results are in Table 1 of Kogut et al. (1996). In detail, 



yn 



the Haslam map was scaled using a power law index a s 
The SFD free-free map was scaled using as = -2.15 and the 
SFD dust model was scaled using an index adust = 2.0. This 
model is known to over-predict the plane emission at DMR 
resolution, thus it is likely to be conservative. In the foreground 
maps, we did not include point sources. Sokasian, Gawiser & 
Smoot (2001) have compiled and analyzed the available extra- 
galactic point source data and have concluded that such sources 



Fig. 8. The histograms of the 94 GHz foreground emission with 
WMAP resolution. Solid, dotted and dashed lines correspond 
to Galactic cut b cut = 5, 10, 20 degrees, respectively. 



Fig. (jSJl shows the histogram of foreground contributions to 
the WMAP 94GHz channel for 3 values of Galactic cut: the 
solid, dotted and dashed lines correspond to \b\ cut = 5°, 10° 
and 20°, respectively. For comparison at v = 2 the value CMB 
contribution to the remaining pixels will be ~ 200^K. Clearly, 
foregrounds are not likely to significantly affect the method. 

Pixel Noise Distribution 



Fig. 9. Histograms of pixel noise variance for the WMAP 94 
GHz channel with pixels size of 7 arcmins. 



4.2.2. Noise model 

Our statistic £ v has been worked out for homogeneous Gaussian 
fields. The inversion of £ v (C(6)/Cq) into the correlation func- 
tion C(6) is accurate provided that signal and noise correla- 
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tion functions are uniform across the sky. For most experi- 
ments this is not the case as the sky coverage is not homo- 
geneous and the noise variance changes with location. For 
this more realistic case we adopt the following strategy: if 
the noise variance (NN) changes according to the number 



of times a pixel has been observed ((TV/TV,) 



p -V2 

l observ 



), pixel 



i will be selected as peak above a threshold v if it verifies 
\6Ti\lv(C CMB (0) + {NiNi)) 112 , i.e., we take into account the 
local noise variation. In the expression, C CMB (0) is the vari- 
ance of the cosmological signal. We apply this formalism to 
the noise model of the WMAP two-year scanning strategy. This 
model gives the number of observations (N \, s ) in each pixel at 
the end of the second mission year. Fig. (|9jl shows the distribu- 
tion of the r.m.s. noise distribution for WMAP 94 GHz chan- 
nel, for 7 arcmin-sized pixels. In the figure, the noise average 
is ^{NN) fu n S ky - per pixel. 

Fig. fllUi shows that this gives similar results as for homo- 
geneous noise. The thick solid line corresponds to the theoreti- 
cal estimate of £ y= 2 according to eq l|3}. For a fixed cosmolog- 
ical signal, if we add homogeneous noise to a CMB map, the 
dot-dashed line will be the numerical estimate, like in Fig. Q. 
The thin solid line gives the estimate of £ v= 2 assuming the noise 
is homogeneous when it is not, while the long-dashed line gives 
the estimate when the pixels have been selected according to 
the criterion defined in the previous paragraph. The latter result 
is comparable to the case of homogeneous noise and, in this 
sense, this procedure gives an optimal estimator. 



4.2.3. Results 

We performed 150 simulations of WMAP data using the 
SCDM model, with inhomogeneous noise and foreground 
residuals. We choose three different galactic cuts and two dif- 
ferent amplitudes for the foreground residuals: 10% and 20% 
amplitude of the original foreground data at 94GHz. In the first 
case, we performed simulations with a Galactic cut at \b\ = 5° 
and 10° and for the 20% residuals we imposed a cut at \b\ — 20°. 
We have checked that these models give a foreground residual 
level comparable to the actual contribution found in the data by 
Bennett et al. (2003c). 




Fig. 11. Radiation power spectra obtained from maps contain- 
ing inhomogeneous noise and foreground residuals. The cor- 
relation function was computed using all pixels above v = 2. 
As in Fig. @, solid lines correspond to the input power spec- 
trum. The dotted lines give the average of 150 simulations. The 
dashed lines show the l-cr confidence level given by the op- 
timal variance of Eq. The shaded area corresponds to the 
l-cr error region obtained from the simulations. 




8 (degrees) 

Fig. 10. Thick solid and dot-dashed lines: £ v= 2 of the input 
model and estimator when the noise is homogeneous, as in 
Fig. ([0. When the noise is inhomogeneous, the thin solid line 
gives % v= 2 obtained when we apply the same method as if the 
noise were homogeneous. The long-dashed line corresponds to 
the case when we take into account the variable amplitude of 
the noise and select pixels accordingly (see text). The shaded 
area is the 1-sigma optimal variance error bars. 



In Fig. ( II 11 we show the results of 150 simulations. As 
in Fig. ®, the dotted lines shows the mean power spectrum 
of all the simulations, while the solid line shows the input 
model. Shaded areas give the l-cr confidence level and dashed 
lines the optimal variance l-cr error bars (see Eq. (jfji). In the 
whole range, the error bars were almost identical (less than 5% 
increase) to the case with no foregrounds and homogeneous 
noise. The correlations in /-space are practically the same as 
before (Al 10). The deviations of the mean from the input 
model are from the smaller number of simulations. 

To quantify the amplitude of the non-Gaussian compo- 
nents, we computed the skewness and kurtosis of our simu- 
lations. The skewness and kurtosis are defined (cf. Press et 
al. 1992) as cr 3 = ((6T/o- 2 )\ cr A = ((6T/o- 2 ) 4 } - 3, respec- 
tively, with o"2 being the variance of the data. In our case, the 
non-gaussian signal is dominated by the foreground residuals. 
We computed the skewness and kurtosis for a randomly se- 
lected map of the previous set of simulations. For a Galactic 
cut at 10° and 10% foreground residuals amplitude, the values 
found were cr 3 = 0.12 and cr 4 = 1.2. For comparison, for the 
same CMB map re-simulated with a realistic noise component, 
galactic cut and no foregrounds, the values were cr 3 = 0.02 and 
cr 4 = -0.03. 

To conclude, the addition of a non-Gaussian signal (the 
foreground residuals) and inhomogeneous noise does not affect 
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significantly the performance of the method. The reason that 
the introduction of Galactic foregrounds produces only small 
variations in the results is at the very core of our procedure: 
with the new KHA statistic, Eq.[3]we select regions with high 
S/N ratio, the peaks, and small departures from gaussianity nei- 
ther degrades the quality of our subset of data, nor introduces 
additional correlations. 

The results presented in Fig. dl 11 1 would improve for lower 
thresholds. As the method provides correlated C/'s, it is nec- 
essary to bin the estimates into band powers. Other methods 
share, to some extent, this limitation; direct computation of the 
power spectrum from a map which necessarily has a cut due to 
bright foreground regions, is possible only with a finite band- 
width determined by size of the Galactic cut. This was avoided 
in the COBE/DMR data with the Gramm-Schmidt orthogonal- 
ization of the base functions (Gorski 1994), but the procedure is 
impractical for the mega-pixel CMB datasets. For the WMAP 
94 GHz channel, our method provides, under reasonable ac- 
curacy, around 12-14 independent data points for the first two 
Doppler peaks of the CMB power spectrum, and with a gain 
factor of 10 3 - 10 4 in CPU time compared to standard meth- 
ods. 

We have also demonstrated that our method retains its 
accuracy in the presence of foregrounds and for the realis- 
tic/inhomogeneous WMAP noise. The treatment in the pres- 
ence of the inhomogeneous noise can also be improved by 
removing the pixels with significantly fewer observations. 
Furthermore, because the peaks method computes a correla- 
tion function, £ v , it is immune to masking. The latter would 
allow us to reduce the foregrounds signal by removing from 
the CMB maps more isolated regions with higher foreground 
contribution. 



5. Application to Planck 

The Planck satellite 2 is due to be launched in early 2007 to map 
the all sky CMB distribution with an even finer resolution than 
WMAP. It will have two instruments: the low-frequency instru- 
ment (LFI), that operates at three frequency channels between 
30 and 70 GHz, and the high-frequency instrument (HFI) in 
six frequency channels between 100 and 857 GHz. The high- 
est Planck resolution will be 5 arcmin and the instrument noise 
after one year of operations is expected to be around 6-10 //K. 

For this case we simulated the 217 GHz channel, with a 
beam of 5.5' FWHM and an average noise level of 1 1.4//K per 
beam area. We tested this channel in two different cosmological 
models: for consistency, we performed 125 simulations with 
the same SCDM model used for WMAP, (top row of Fig. ( fT2l . 
but we also tested our method with the ACDM concordance 
model, given by Q m = 0.292, Q A = 0.708, Q b = 0.044 and 
h = 0.72, (bottom row of the same figure). We assumed that 
the instrumental noise was the sum of two different components 
with white noise and 1// contributions (see Maino et al. 1999 
for a detailed account of systematic effect on the Planck LFI 
instrument). 



2 I http://astro.estec.esa.nl/Planck/ 1 



The 1 // component is characterized by a knee frequency fa, 
for which the power spectra of both white and 1 // noise com- 
ponents are equal. If the spin frequency fa is not much smaller 
than the knee frequency fa, then we can use the fact that the 
telescope spends sixty minutes in each ring to simply consider 
the average noise pattern for every ring, (Maino et al. 2002). In 
other words, we can neglect the noise high frequency compo- 
nents present in the same ring, as these will be averaged out. 
We model the 1 // component as a different baseline present 
in every scanned ring, giving a pattern of stripes in the noise 
map (Janssen et al. 1996). We used the Planck 217 GHz chan- 
nel noise model, provided to us by the Max Planck Institute 
fur Astrophysik (MPA) at Garching, that includes all system- 
atic effects. For each simulation, we added a realization of this 
1 // component to the cosmological and white noise compo- 
nents. For simplicity we did not include foreground residuals 
as in the previous section we have shown to give negligible 
contribution. For similar reasons we also did not include point 
sources. 

The maximum multipole to which a beam of 5.5' FWHM 
is sensitive is /b ea m = 1472. The HEALPix configuration cho- 
sen to pixelize the CMB sky as seen by this channel was 
Miide = 1024, which yields pixels of 3.5 arcmins. In Fig. d!2i 
we show the results of 500 simulations. As indicated, we have 
chosen three threshold levels that have the same number of pix- 
els that the three v's considered for WMAP. We plot the results 
for those thresholds following the conventions of Fig. Note 
that, while we obtain the radiation power spectrum from the 
same number of pixels as in the case of WMAP, the method 
recovers it up to / = 1472, the largest possible multipole de- 
termined by the beam size. This could be expected, since the 
clustering of peaks and troughs is a statistic particularly sensi- 
tive to small angular scales. In the range e [#beam, ~ 10°], or 
equivalently / e [40, /beam], our method reconstructs the power 
spectrum with an accuracy depending only on the threshold. 
The range in /-space probed by the method depends on the res- 
olution of each experiment but is independent of the total num- 
ber of peaks used in the analysis. 

Using the data from the three set of simulations we can 
compute how the band power error bars scale with the threshold 
v. In this case, a fit of the form <x C( oc iVf eaks to the error bars of 
all multipoles gives: (3 = -0.65 + 0.011, close to the expected 
behavior cr Cl oc 1 / -y//V pea ks ■ As mentioned before, this scaling 
can be used to decide what threshold to choose to achieve a 
preselected accuracy. The correlations on the C/'s were similar 
to those of WMAP, of the order of Al ^ 10. 

6. Discussion and conclusions 

We have presented the detailed implementation of the peaks 
KHA method for computing the CMB power spectrum. The 
method is based on the correlation properties of peaks on the 
CMB and assumes that temperature fluctuations are Gaussian. 
We have shown, that the method is robust in the presence of 
Galactic (non-Gaussian) foregrounds and realistic WMAP and 
Planck noise inhomogeneities. 

This method provides a new and independent way to char- 
acterize the sub-degree structure of the CMB. Its accuracy is 
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Fig. 12. Radiation power spectrum for the Planck 217 GHz 
channel: in the top row we have considered the SCDM model 
used in WMAP simulations. In the top panels we follow the 
conventions of Fig. ®: solid lines represent the input model, 
dotted lines the mean of 125 simulations, dashed lines are opti- 
mal variance error bars at lcr and the shaded areas are the same 
confidence limit obtained from the simulations. In the bottom 
row, we show results for the concordance ACDM model for 
one single realization (dotted line). Again, the input model is 
given by the solid line and the optimal variance error bars are 
given, as before, by the dashed lines. 
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Fig. 13. The value of the peak threshold where the number 
of operations in the KHA method equals AT 1 ' 61 , is plotted as a 
function of the number of pixels. 



generally lower than standard methods, but it is remarkably 
faster for the range of values of v required by WMAP and 
Planck. This is demonstrated in Fig. dl3> . where we compare 
speed of a standard 0(N log N) ~ N 161 method with ours. 

For WMAP the KHA method should work most optimally 
for peak thresholds of v ~ 1.8 - 2.5 or in only (2.5 x 10 3 - 
8 x 10~ 5 )N 2 operations. Applying the method to the simulated 




i/=1.6, 6x10 3 N pi™] operations 

..I i i I.. 

300 400 500 600 
Multipole / 



Fig. 14. Radiation power spectrum that would be recovered us- 
ing the WMAP second year data with a threshold of v = 1.6: 
it was obtained in only 6 x 10~ 3 N 2 operations. The solid line 
corresponds to the input model and the dotted line is the raw 
power spectrum recovered from one single simulation, prior to 
binning. The shaded area shows the extrapolated l-cr uncer- 
tainties of the CMB power spectrum binned with Al — 45. 



PLANCK sky maps has shown that we could probe, with the 
same number of points, the power spectrum out to much higher 
multipoles and without significant loss of accuracy. For Planck 
mission parameters, the KHA method should work most op- 
timally for v ~ 2 - 3, reducing the number of steps down to 
~ 4 x 10 _6 A^ 2 operations for thresholds as high as v ~ 3. This 
can potentially lead to a gain of up to ~ 10 5 in speed compared 
to the existing 0(N 2 ) methods. 

A significant advantage is that because we use a £ v statistic, 
the computed CMB spectrum is immune to masking, allowing 
to remove high noise and foreground emission parts of the sky. 
This also makes straightforward the analysis of particular iso- 
lated regions in the sky, enabling its application to future CMB 
experiments of very high angular resolution scanning small ce- 
lestial patches. 

With the scalings of <x C( with v for both WMAP and Planck 
one can estimate the uncertainties of the different multipoles 
at still lower peak threshold values. In Fig. (I14> we show what 
would be the spectrum recovered from the 2 year WMAP data 
release using a threshold v - 1.6. The whole calculation would 
require only 6 X lO -3 ^ 2 operations. 
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